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Abstract 

Inpainting,  the  technique  of  modifying  an  image  in  an  undetectable 
form,  is  as  ancient  as  art  itself.  The  goals  and  applications  of  in¬ 
painting  are  numerous,  from  the  restoration  of  damaged  paintings 
and  photographs  to  the  removal/replacement  of  selected  objects.  In 
this  paper,  we  introduce  a  novel  algorithm  for  digital  inpainting  of 
still  images  that  attempts  to  replicate  the  basic  techniques  used  by 
professional  restorators.  After  the  user  selects  the  regions  to  be 
restored,  the  algorithm  automatically  fills-in  these  regions  with  in¬ 
formation  surrounding  them.  The  fill-in  is  done  in  such  a  way  that 
isophote  lines  arriving  at  the  regions  boundaries  are  completed  in¬ 
side.  In  contrast  with  previous  approaches,  the  technique  here  in¬ 
troduced  does  not  require  the  user  to  specify  where  the  novel  in¬ 
formation  comes  from.  This  is  automatically  done  (and  in  a  fast 
way),  thereby  allowing  to  simultaneously  fill-in  numerous  regions 
containing  completely  different  structures  and  surrounding  back¬ 
grounds.  In  addition,  no  limitations  are  imposed  on  the  topology  of 
the  region  to  be  inpainted.  Applications  of  this  technique  include 
the  restoration  of  old  photographs  and  damaged  film;  removal  of  su¬ 
perimposed  text  like  dates,  subtitles,  or  publicity;  and  the  removal 
of  entire  objects  from  the  image  like  microphones  or  wires  in  spe¬ 
cial  effects. 
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1  Introduction 

The  modification  of  images  in  a  way  that  is  non-detectable  for  an 
observer  who  does  not  know  the  original  image  is  a  practice  as  old 
as  artistic  creation  itself.  Medieval  artwork  started  to  be  restored  as 
early  as  the  Renaissance,  the  motives  being  often  as  much  to  bring 
medieval  pictures  “up  to  date”  as  to  fill  in  any  gaps  [1,  2].  This 
practice  is  called  retouching  or  inpainting.  The  object  of  inpainting 
is  to  reconstitute  the  missing  or  damaged  portions  of  the  work,  in 
order  to  make  it  more  legible  and  to  restore  its  unity  [2]. 

The  need  to  retouch  the  image  in  an  unobtrusive  way  extended 
naturally  from  paintings  to  photography  and  film.  The  purposes 
remain  the  same:  to  revert  deterioration  (e.g.,  cracks  in  photographs 
or  scratches  and  dust  spots  in  film),  or  to  add  or  remove  elements 
(e.g.,  removal  of  stamped  date  and  red-eye  from  photographs,  the 
infamous  “airbrushing”  of  political  enemies  [3]). 

Digital  techniques  are  starting  to  be  a  widespread  way  of  inpaint¬ 
ing,  ranging  from  attempts  to  fully  automatic  detection  and  removal 
of  scratches  in  film  [4,  5,  6],  all  the  way  to  software  tools  that  allow 
a  sophisticated  but  mostly  manual  process  [7]. 

In  this  article  we  introduce  a  novel  algorithm  for  automatic  digi¬ 
tal  inpainting,  being  its  main  motivation  to  replicate  the  basic  tech- 
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niques  used  by  professional  restorators.  At  this  point,  the  only  user 
interaction  required  by  the  algorithm  here  introduced  is  to  mark 
the  regions  to  be  inpainted.  Although  a  number  of  techniques  ex¬ 
ist  for  the  semi-automatic  detection  of  image  defects  (mainly  in 
films),  addressing  this  is  out  of  the  scope  of  this  paper.  Moreover, 
since  the  inpainting  algorithm  here  presented  can  be  used  not  just 
to  restore  damaged  photographs  but  also  to  remove  undesired  ob¬ 
jects  and  writings  on  the  image,  the  regions  to  be  inpainted  must  be 
marked  by  the  user,  since  they  depend  on  his/her  subjective  selec¬ 
tion.  Here  we  are  concerned  on  how  to  “fill-in”  the  regions  to  be 
inpainted,  once  they  have  been  selected.1  Marked  regions  are  au¬ 
tomatically  filled  with  the  structure  of  their  surrounding,  in  a  form 
that  will  be  explained  later  in  this  paper. 

2  Related  work  and  our  contribution 

We  should  first  note  that  classical  image  denoising  algorithms  do 
not  apply  to  image  inpainting,  since  the  regions  to  be  inpainted  are 
usually  large.  That  is,  regions  occupied  by  top  to  bottom  scratches 
along  several  film  frames,  long  cracks  in  photographs,  superim¬ 
posed  large  fonts,  and  so  on,  are  of  significant  larger  size  than  the 
type  of  noise  assumed  in  common  image  enhancement  algorithms. 
In  addition,  in  common  image  enhancement  applications,  the  pixels 
contain  both  information  about  the  real  data  and  the  noise  (e.g.,  im¬ 
age  plus  noise  for  additive  noise),  while  in  image  inpainting,  there 
is  no  significant  information  in  the  region  to  be  inpainted.  The 
information  is  mainly  in  the  regions  surrounding  the  areas  to  be 
inpainted.  There  is  then  a  need  to  develop  specific  techniques  to 
address  these  problems. 

Mainly  three  groups  of  works  can  be  found  in  the  literature  re¬ 
lated  to  digital  inpainting.  The  first  one  deals  with  the  restoration 
of  films,  the  second  one  is  related  to  texture  synthesis,  and  the  third 
one,  a  significantly  less  studied  class  though  very  influential  to  the 
work  here  presented,  is  related  to  disocclusion. 

Joyeux  et  al.  [4]  devised  a  2-steps  frequency  based  reconstruc¬ 
tion  system  for  detecting  and  removing  line  scratches  in  films.  They 
propose  to  first  recover  low  and  then  high  frequencies.  Although 
good  results  are  obtained  for  their  specific  application,  the  algo¬ 
rithm  can  not  handle  large  loss  areas.  Frequency  domain  algorithms 
trade  a  good  recovery  of  the  overall  structure  of  the  regions  for  poor 
spatial  results  regarding,  for  instance,  the  continuity  of  lines. 

Kokaram  et  al.  [6]  use  motion  estimation  and  autoregressive 
models  to  interpolate  losses  in  films  from  adjacent  frames.  The 
basic  idea  is  to  copy  into  the  gap  the  right  pixels  from  neighboring 
frames.  The  technique  can  not  be  applied  to  still  images  or  to  films 
where  the  regions  to  be  inpainted  span  many  frames. 

Hirani  and  Totsuka  [8]  combine  frequency  and  spatial  domain  in¬ 
formation  in  order  to  fill  a  given  region  with  a  selected  texture.  This 
is  a  very  simple  technique  that  produces  incredible  good  results.  On 
the  other  hand,  the  algorithm  mainly  deals  with  texture  synthesis 


lln  order  to  study  the  robustness  of  the  algorithm  here  proposed,  and  not 
to  be  too  dependent  on  the  marking  of  the  regions  to  be  inpainted,  we  mark 
them  in  a  very  rough  form  with  any  available  paintbrush  software.  Marking 
these  regions  in  the  examples  reported  in  this  paper  just  takes  a  few  seconds 
to  a  non-expert  user. 
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(and  not  with  structured  background),  and  requires  the  user  to  select 
the  texture  to  be  copied  into  the  region  to  be  inpainted.  For  images 
where  the  region  to  be  replaced  covers  several  different  structures, 
the  user  would  need  to  go  through  the  tremendous  work  of  segment¬ 
ing  them  and  searching  corresponding  replacements  throughout  the 
picture.  Although  part  of  this  search  can  be  done  automatically, 
this  is  extremely  time  consuming  and  requires  the  non-trivial  selec¬ 
tion  of  many  critical  parameters,  e.g.,  [9].  Other  texture  synthesis 
algorithms,  e.g.,  [9,  10,  11],  can  be  used  as  well  to  re-create  a  pre¬ 
selected  texture  to  fill-in  a  (square)  region  to  be  inpainted. 

In  the  group  of  disocclusion  algorithms,  a  pioneering  work  is 
described  in  [12].  The  authors  presented  a  technique  for  removing 
occlusions  with  the  goal  of  image  segmentation.2  The  basic  idea 
is  to  connect  T-junctions  at  the  same  gray-level  with  elastica  min¬ 
imizing  curves.  The  technique  was  mainly  developed  for  simple 
images,  with  only  a  few  objects  with  constant  gray-levels,  and  will 
not  be  applicable  for  the  examples  with  natural  images  presented 
later  in  this  paper.  Masnou  and  Morel  [13]  recently  extended  these 
ideas,  presenting  a  very  inspiring  general  variational  formulation 
for  disocclusion  and  a  particular  practical  algorithm  implementing 
some  of  the  ideas  in  this  formulation.  The  algorithm  performs  in¬ 
painting  by  joining  with  geodesic  curves  the  points  of  the  isophotes 
arriving  at  the  boundary  of  the  region  to  be  inpainted.  As  reported 
by  the  authors,  the  regions  to  be  inpainted  are  limited  to  having 
simple  topology,  e.g.,  holes  are  not  allowed.3  In  addition,  the  angle 
with  which  the  level  lines  arrive  at  the  boundary  of  the  inpainted 
region  is  not  (well)  preserved:  the  algorithm  uses  straight  lines  to 
join  equal  gray  value  pixels.  These  drawbacks  will  be  exemplified 
later  in  this  paper.  On  the  other  hand,  we  should  note  that  this  is  the 
closest  technique  to  ours  and  has  motivated  in  part  and  inspired  our 
work. 

2.1  Our  contribution 

Algorithms  devised  for  film  restoration  are  not  appropriate  for  our 
application  since  they  normally  work  on  relatively  small  regions 
and  rely  on  the  existence  of  information  from  several  frames. 

On  the  other  hand,  algorithms  based  on  texture  synthesis  can  fill 
large  regions,  but  require  the  user  to  specify  what  texture  to  put 
where.  This  is  a  significant  limitation  of  these  approaches,  as  may 
be  seen  in  examples  presented  later  in  this  paper,  where  the  region 
to  be  inpainted  is  surrounded  by  hundreds  of  different  backgrounds, 
some  of  them  being  structure  and  not  texture. 

The  technique  we  propose  does  not  require  any  user  interven¬ 
tion,  once  the  region  to  be  inpainted  has  been  selected.  The  algo¬ 
rithm  is  able  to  simultaneously  fill  regions  surrounded  by  different 
backgrounds,  without  the  user  specifying  “what  to  put  where.”  No 
assumptions  on  the  topology  of  the  region  to  be  inpainted,  or  on 
the  simplicity  of  the  image,  are  made.  The  algorithm  is  devised 
for  inpainting  in  structured  regions  (e.g.,  regions  crossing  through 
boundaries),  though  it  is  not  devised  to  reproduce  large  textured 
areas.  As  we  will  discuss  later,  the  combination  of  our  proposed 
approach  with  texture  synthesis  techniques  is  the  subject  of  current 
research. 

3  The  digital  inpainting  algorithm 

3.1  Fundamentals 

Let  Q  stand  for  the  region  to  be  inpainted,  and  <90  for  its  boundary 
(note  once  again  that  no  assumption  on  the  topology  of  O  is  made). 

2 Since  the  region  to  be  inpainted  can  be  considered  as  occluding  objects, 
removing  occlusions  is  analogous  to  image  inpainting. 

3  This  is  not  intrinsic  to  the  general  variational  formulation  they  propose, 

only  to  the  specific  discrete  implementation  they  perform. 


Figure  1:  ID  signals  Io(i)  and  b(i). 


Intuitively,  the  technique  we  propose  will  prolong  the  isophote  lines 
arriving  at  <90,  while  maintaining  the  angle  of  “arrival.”  We  pro¬ 
ceed  drawing  from  <90  inward  in  this  way,  while  curving  the  pro¬ 
longation  lines  progressively  to  prevent  them  from  crossing  each 
other.4 

Before  presenting  the  detailed  description  of  this  technique,  let 
us  analyze  how  experts  inpaint.  Conservators  of  the  Minneapolis 
Institute  of  Arts  were  consulted  for  this  work  and  made  it  clear  to 
us  that  inpainting  is  a  very  subjective  procedure,  different  for  each 
work  of  art  and  for  each  professional.  There  is  no  such  thing  as 
“the”  way  to  solve  the  problem,  but  the  underlying  methodology  is 
as  follows:  (1.)  The  global  picture  determines  how  to  fill  in  the  gap, 
the  purpose  of  inpainting  being  to  restore  the  unity  of  the  work;  (2.) 
The  structure  of  the  area  surrounding  O  is  continued  into  the  gap, 
contour  lines  are  drawn  via  the  prolongation  of  those  arriving  at 
dO;  (3.)  The  different  regions  inside  O,  as  defined  by  the  contour 
lines,  are  filled  with  color,  matching  those  of  dO;  and  (4.)  The  small 
details  are  painted  (e.g.  little  white  spots  on  an  otherwise  uniformly 
blue  sky).  In  other  words,  “texture”  is  added. 

Our  algorithm  simultaneously  and  iteratively  performs  the  steps 
(2.)  and  (3.)  above.  (In  the  discussion  section  we  will  argue  how 
both  steps  can  be  performed  separately,  and  we  will  also  discuss 
step  (4.).)  We  progressively  “shrink”  the  gap  O  by  prolonging  in¬ 
ward,  in  a  smooth  way,  the  lines  arriving  at  the  gap  boundary  <90. 

3.2  A  one  dimensional  illustration 

A  simple  one  dimensional  (ID)  example  is  first  used  to  illustrate 
our  technique.  This  is  done  only  for  a  didactic  purpose,  since  our 
algorithm  was  devised  for  2D,  and  there  are  other  techniques  (such 
as  splines)  that  might  yield  better  interpolation  results  in  ID.  Our 
ID  example  will  help  to  explain  how  we  translate  the  basic  inpaint¬ 
ing  ideas  presented  above  into  a  working  algorithm.  Consider  a  ID 
(discrete)  signal  Io(i)  :  [0,  M]  — >  M  (  [0,  M\  C  IV  )  and  a  “box- 
function”  mask  signal  b(i),  as  illustrated  in  Figure  1.  The  region 
O  to  inpaint  is  given  by  the  set  of  points  where  b(i)  ^  0,  while 
its  boundary  dO  are  the  points  Q  and  P.  Let  be  the  direction  of 
increasing  i. 

The  progressive  nature  of  the  inpainting  procedure  suggests  an 
“evolution”  approach.  We  define  a  family  of  ID  signals  I(i,  n)  : 
[0,  M\  x  [0,  oo)  — >  M  ( [0,  M\  C  IV,  [0,  oo]  C  IV )  where  n  is  the 
“time”  parameter,  such  that  /(i,  0)  =  Io(i)  and  Umn^ooI(i,  n )  — 
Ir(i),  where  Ir  (i)  is  the  output  of  the  algorithm  (inpainted  image). 
At  each  iteration  step  n  we  want  to  construct  I(i,n)  given  I(i,n  — 
!)• 

For  the  ID  case,  the  “prolongation  of  lines  into  O”  is  interpreted 
as  “the  continuation  of  the  function  Io(i)  from  the  left  of  Q  to  its 
right,  and  from  the  right  of  P  to  its  left,  in  a  smooth  way.”  Therefore, 


4This  type  of  information  propagation  is  related  and  might  be  applicable 
to  velocity  fields  extension  in  level-set  techniques  [14]. 


Figure  2:  Discrete  example;  see  text.. 

in  order  to  inpaint  in  the  direction  of  7$,  the  following  pseudo-code 
realizes  the  ideas  above: 

for  n  =  0  to  n  =  Number  —  of  —  Iterations 
for  each  point  (pixel)  (i)  E  Q 

compute  the  smoothness  of  I71  ( i ) :  Ln  ( i )  :=  Ixx  (i) 

compute  the  smoothness  of  In  (i  —  1):  Ln  ( i  —  1)  :=  Ixx  ( i  —  1) 

compute  the  smoothness  variation  in  —  i\h 

1) 

r+1(i)  =  r(i)  +  Ati?(i) 
end 
end, 

where  the  index  i  denotes  spatial  position,  the  superindex  n  de¬ 
notes  the  artificial  evolution  time  ( I(i ,  n )  =  the  subindex 

x  denotes  (discretized)  spatial  differentiation  in  the  horizontal  axis 
direction,  and  At  controls  the  velocity  of  the  evolution.  We  are 
changing  in  time  the  function  I(i,n)  so  that  at  each  point  i,  In+1  ( i ) 
tends  to  follow  In  (t—  1) .  The  purpose  of  the  smoothness  estimation 
is  to  ensure  continuity  not  only  of  I (i,  n)  but  also  of  its  derivatives. 
Our  discrete  smoothness  estimator  is  the  second  spatial  derivative 
of  I(i ,  n),  Ln(i )  =  IxX(i ),  which  in  discrete  form  can  be  written 
as  Ln(i )  =  T(i-  1)  -  2 In{i)  +In(i  +  1). 

To  explain  why  this  works,  consider  Figure  2.  The  continuous 
line  denotes  the  part  of  I(i,  n)  that  is  not  changed  by  the  algorithm 
since  it  is  outside  O.  To  the  right  of  Q,  the  three  dashed  lines  denote 
three  possible  initial  conditions  for  I(x,  0)  inside  O.  Being  straight 
lines,  L(R ,  0)  =  0  in  the  three  cases.  At  the  point  Q  on  the  other 
hand,  the  value  of  L(Q,  0)  is  different  in  each  case: 

If  I(R,0}  =  a,  L(Q,  0)  >  0,  It(R,0)  =  L{R,  0)  -  L(Q,  0)  <  0,  and 
I(R,  0)  decreases; 

If  I(R,0}  =  b,L(Q,  0)  -  0 0)  =  L(R,  0)  —  L(Q,  0)  -  0,and/(.R,  0) 
does  not  change; 

If  I(R,  0)  =  c,L(Q,  0)  <  0,It(R,0)  =  L(R,0)-L(Q,0)  >  0,andI(R,0) 
increases. 

In  all  three  cases  the  evolution  is  deforming  I(i,  n)  inside  O  into 
a  continuation  of  Io(i)  to  the  left  of  Q. 

Figure  3  shows  three  stages  of  the  evolution  of  the  signal  in  Fig¬ 
ure  1.  The  dotted  line  corresponds  to  the  initial  condition,  I(i ,  0). 
The  dashed  line  corresponds  to  an  intermediate  state.  The  contin¬ 
uous  line  corresponds  to  Ir{i),  the  steady  state  of  the  evolution. 
Notice  how  continuity  up  to  the  first  derivative  is  obtained  for  the 
point  Q.  That  is  not  the  case  for  P,  since  the  direction  of  propaga¬ 
tion  is  7$  and  not  —7$.  We  will  later  see  that  this  is  not  a  problem 
for  the  2D  case.  We  have  tested  ID  examples  with  different  initial 
conditions,  and  found  that  the  algorithm  always  converges  to  the 
same  solution  within  a  0.1%  error  margin.  The  formal  study  of  this 
experimental  result  is  of  theoretical  interest  (see  Section  5). 

3.3  The  2D  inpainting  algorithm 

Let  us  re-write  the  ID  discrete  equation  describing  the  rate  of 
change  of  I(i,  n): 

r+1(o  =  /n(o  +  Ati?(t))  (i) 


Figure  3:  Three  stages  of  the  evolution  of  I(i ,  n)  in  Figure  1. 


Figure  4:  Propagation  direction  as  the  normal  to  the  signed  distance 
to  the  boundary  of  the  region  to  be  inpainted. 


where 

I?(i)  =  Ln(i)-Ln(i-  1),  Vi  GO,  (2) 

Ln(i)  =  In(i  -  1)  -  2 In(i)  +  In(i  +  1).  (3) 

In  other  words,  we  estimate  the  smoothness  Ln(i)  of  our  func¬ 
tion  and  compute  its  change  along  the  —  1$  direction,  adapting  the 
image  in  the  region  to  be  inpainted  with  this  variation  of  smooth¬ 
ness  quantity.  Let  us  now  see  how  to  extend  these  ideas  into  2D. 

Let  Io(i,j)  :  [0,  M]  x  [0,  N]  JR,  with  [0,  M]  x  [0,  N]  C 
IV  x  IN,  be  a  discrete  2D  gray  level  image.  The  inpainting  pro¬ 
cedure  will  construct  a  family  of  images  I(i,j,n )  :  [0,  M\  x 
[0,iV]  x  [0,  oo)  — >  M  such  that  0)  =  and 

Umn^ooI(ii  j,  n)  —  /j?(i,j),  where  is  the  output  of  the 

algorithm  (inpainted  image). 

Following  the  ID  procedure,  first  we  estimate  the  smoothness 
Ln(i,j)  of  our  function.  For  this  purpose  we  may  use  a  simple 
discrete  implementation  of  the  Laplacian:  Ln(i,j)  :=  I£x(h  j)  + 
Iyy(i ,  j).  Other  smoothness  estimators  might  be  used,  though  sat¬ 
isfactory  results  were  already  obtained  with  this  very  simple  selec¬ 
tion. 

Then,  we  must  compute  the  change  of  this  value  along  —7$.  In 
order  to  do  this  we  must  first  define  what  the  direction  7$  for  the 
2D  information  propagation  will  be.  One  possibility  is  to  define  N 
as  the  normal  to  the  signed  distance  to  dfl,  i.e.,  at  each  point  (i,  j) 
in  O  the  vector  7$ (i,j)  will  be  normal  to  the  “shrinked  version”  of 
d£l  to  which  (i,  j)  belongs,  see  Figure  4.  This  choice  is  motivated 
by  the  belief  that  a  propagation  normal  to  the  boundary  would  lead 
to  the  continuity  of  the  isophotes  at  the  boundary.  Instead,  what 
happens  is  that  the  lines  arriving  at  dfl  curve  in  order  to  align  with 
7$,  see  Figure  5.  This  is  of  course  not  what  we  expect.  Note  that 
the  orientation  of  dfl  is  not  intrinsic  to  the  image  geometry,  since 
the  region  to  be  inpainted  !  !  is  arbitrary. 

If  isophotes  tend  to  align  with  7$,  the  best  choice  for  7$  is  then 
the  isophotes  directions.  This  is  a  bootstrapping  problem:  hav¬ 
ing  the  isophotes  directions  inside  O  is  equivalent  to  having  the 
inpainted  image  itself,  since  we  can  easily  recover  the  gray  level 
image  from  its  isophote  direction  field  (  see  the  discussion  section 
and  [15]). 


Figure  5:  Unsuccessful  choice  of  the  information  propagation  di¬ 
rection.  Left:  detail  of  the  original  image,  region  to  be  inpainted  is 
in  white.  Right:  restoration. 


describe.  Every  few  iterations,  a  step  of  anisotropic  diffusion  is  ap¬ 
plied  (a  straightforward,  central  differences  implementation  of  (4) 
is  used;  for  details  see  [16,  17]).  This  process  is  repeated  until  a 
steady  state  is  achieved. 

Let  In(i,j)  stand  for  each  one  of  the  image  pixels  inside  the 
region  O  at  the  inpainting  “time”  n.  Then,  the  discrete  inpainting 
equation  borrows  from  the  numerical  analysis  literature  and  is  given 
by 

where 


We  use  then  a  time  varying  estimation  of  the  isophotes  direction 
field:  for  any  given  point  (i,  j),  the  (discretized)  gradient  vector 
VIn(i,j)  gives  the  direction  of  largest  (spatial)  change,  while  its 
90  degrees  rotationV±/n(i,  j)  is  the  direction  of  smallest  (spatial) 
change,  so  the  vector  S/±In  ( i ,  j)  gives  the  isophotes  direction.  Our 
field  A^  is  then  given  by  the  time- varying  A^(i,  j,  n)  =  S/±In(i,  j). 
We  are  using  a  time- varying  estimation  that  is  coarse  at  the  begin¬ 
ning  but  progressively  achieves  the  desired  continuity  at  dO,  in¬ 
stead  of  a  fixed  field  7$(i,j)  that  would  imply  to  know  the  direc¬ 
tions  of  the  isophotes  from  the  start. 

Notethat  the  field  is  not  normalized  as  it  was  in  the  ID  case, 
where  N  was  just  the  direction  of  increasing  i.  The  reason  for  this 
choice  relies  on  the  numerical  stability  of  the  algorithm,  and  will  be 
discussed  in  the  following  subsection. 

Since  we  are  performing  an  inpainting  along  the  isophotes,  it  is 
irrelevant  if  S7±In  ( i ,  j)  is  obtained  as  a  clockwise  or  counterclock¬ 
wise  rotation  of  VIn(i ,  j).  In  both  cases,  the  change  of  In(i ,  j ) 
along  those  directions  should  be  minimum. 

Recapping,  we  estimate  a  variation  of  the  smoothness,  given  by 
a  discretization  of  the  2D  Laplacian  in  our  case,  and  project  this 
variation  into  the  isophotes  direction.  This  projection  is  used  to 
update  the  value  of  the  image  inside  the  region  to  be  inpainted. 

To  ensure  a  correct  evolution  of  the  direction  field,  a  diffusion 
process  is  interleaved  with  the  image  inpainting  process  described 
above.  That  is,  every  few  steps  (see  below),  we  apply  a  few  itera¬ 
tions  of  image  diffusion.  This  diffusion  corresponds  to  the  periodi¬ 
cal  curving  of  lines  to  avoid  them  from  crossing  each  other,  as  was 
mentioned  in  section  3.1.  We  use  anisotropic  diffusion,  [16,  17], 
in  order  to  achieve  this  goal  without  losing  sharpness  in  the  recon¬ 
struction.  In  particular,  we  apply  a  straightforward  discretization  of 
the  following  continuous-time/continuous- space  anisotropic  diffu¬ 
sion  equation: 

-  9c(x,y)n(x,y,t) \VI(x,y,t)\ ,V(x,y)  €  (4) 

where  De  is  a  dilation  of  O  with  a  ball  of  radius  e,  n  is  the  Euclidean 
curvature  of  the  isophotes  of  u  and  gc  ( x ,  y )  is  a  smooth  function  in 
Oe  such  that  g^(x,y)  —  0  in  Oe  \  O,  and  gc(x,y)  —  1  at  the  set  of 
points  of  D  whose  distance  to  dQ  is  larger  that  e  (this  is  a  way  to 
impose  Dirichlet  boundary  conditions  for  the  equation  (4)). 


3.4  Discrete  scheme  and  implementation  details 

The  only  input  to  our  algorithm  are  the  image  to  be  restored  and 
the  mask  that  delimits  the  portion  to  be  inpainted.  As  a  preprocess¬ 
ing  step,  the  whole  original  image  undergoes  anisotropic  diffusion 
smoothing.  The  purpose  is  to  minimize  the  influence  of  noise  on 
the  estimation  of  the  direction  of  the  isophotes  arriving  at  dfl.  After 
this,  the  image  enters  the  inpainting  loop,  where  only  the  values  in¬ 
side  O  are  modified.  These  values  change  according  to  the  discrete 
implementation  of  the  inpainting  procedure,  which  we  proceed  to 
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We  first  compute  the  2D  smoothness  estimation  L  in  (8)  and  the 
isophote  direction  ]^/|A^|  in  (9).  Then  in  (10)  we  compute  j3n , 
the  projection  of  si  onto  the  (normalized)  vector  A^,  that  is,  we 
compute  the  change  of  L  along  the  direction  of  aK  Finally,  we 
multiply  /3n  by  a  slope-limited  version  of  the  norm  of  the  gradient 
of  the  image,  |  Vij,  in  (11).  A  central  differences  realization  would 
turn  the  scheme  unstable,  and  that  is  the  reason  for  using  slope- 
limiters.  The  subindexes  b  and  /  denote  backward  and  forward 
differences  respectively,  while  the  subindexes  m  and  M  denote  the 
minimum  or  maximum,  respectively,  between  the  derivative  and 
zero  (we  have  omitted  the  space  coordinates  (i,j)  for  simplicity); 
see  [18]  for  details.  Finally,  let  us  note  that  the  choice  of  V^/ 
instead  of  a  normalized  version  of  it  (as  was  the  case  in  ID)  allows 
for  a  simpler  and  more  stable  numerical  scheme;  see  [19,  20]. 

When  applying  equations  (5)-(l  1)  to  the  pixels  in  the  border  dfl 
of  the  region  O  to  be  inpainted,  known  pixels  from  outside  this 
region  are  used.  That  is,  conceptually,  we  compute  equations  (5)- 
(1 1)  in  the  region  Oe  (an  e  dilation  of  O),  although  update  the  values 
only  inside  O  (that  is,  (5)  is  applied  only  inside  O).  The  information 
in  the  narrow  band  Oe  —  Q  (isophote  directions  and  gray  values)  is 
propagated  inside  D.  Propagation  of  this  information,  both  gray- 
values  and  isophotes  directions,  is  fundamental  for  the  success  of 
the  algorithm. 

In  the  restoration  loop  we  perform  A  steps  of  inpainting  with 
(5),  then  B  steps  of  diffusion  with  (4),  again  A  steps  of  (5),  and 
so  on.  The  total  number  of  steps  is  T.  This  number  may  be  pre- 
established,  or  the  algorithm  may  stop  when  changes  in  the  image 
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Figure  6:  Relation  between  the  (R,G,B)  color  model  and  the  one 
used  in  this  article,  (p,  sine/),  sin'ip). 


are  below  a  given  threshold.  The  values  we  use  are:  A  —  15 ,  B  — 
2,  A t  —  0.1.  The  value  of  T  depends  on  the  size  of  O.  If  O  is 
of  considerable  size,  a  multiresolution  approach  is  used  to  speed-up 
the  process. 

Color  images  are  considered  as  a  set  of  three  images,  and  the 
above  described  technique  is  applied  independently  to  each  one. 
To  avoid  the  appearance  of  spurious  colors,  we  use  a  color  model 
which  is  very  similar  to  the  LUV  model,  with  one  luminance  and 
two  chroma  components.  See  Figure  6. 

4  Results 

The  CPU  time  required  for  inpainting  depends  on  the  size  of  O.  In 
all  the  color  examples  here  presented,  the  inpainting  process  was 
completed  in  less  than  5  minutes  (for  the  three  color  planes),  using 
non-optimized  C++  code  running  on  a  Pentiumll  PC  (128Mb  RAM, 
300MHz)  under  Linux.  All  the  examples  use  images  available  from 
public  databases  over  the  Internet. 

Figure  7  shows,  on  the  left,  a  synthetic  image  with  the  region 
to  inpaint  in  white.  Here  O  is  large  (30  pixels  in  diameter)  and 
contains  a  hole.  The  inpainted  reconstruction  is  shown  on  the  right. 
Notice  that  contours  are  recovered,  joining  points  from  the  inner 
and  outer  boundaries.  Also,  these  reconstructed  contours  follow 
smoothly  the  direction  of  the  isophotes  arriving  at  dfl. 

Figure  8  shows  a  deteriorated  B&W  image  and  its  reconstruc¬ 
tion.  As  in  all  the  examples  in  this  article,  the  user  only  supplied 
the  “mask”  image,  shown  in  Figure  9.  This  image  was  drawn  man¬ 
ually,  using  a  paintbrush-like  program.  The  variables  were  set  to 
the  values  specified  in  the  previous  section,  and  the  number  of  it¬ 
erations  T  was  set  to  3000.  When  multiresolution  is  not  used,  the 
CPU  time  required  by  the  inpainting  procedure  was  approximately 
7  minutes.  With  a  2-level  multiresolution  scheme,  only  2  minutes 
were  needed.  Observe  that  details  in  the  nose  and  right  eye  of  the 
middle  girl  could  not  be  completely  restored.  This  is  in  part  due  to 
the  fact  that  the  mask  covers  most  of  the  relevant  information,  and 
there  is  not  much  to  be  done  without  the  use  of  high  level  prior  in¬ 
formation  (e.g.,  the  fact  that  it  is  an  eye).  These  minor  errors  can  be 
corrected  by  the  manual  procedures  mentioned  in  the  introduction, 
and  still  the  overall  inpainting  time  would  be  reduced  by  orders  of 
magnitude.  In  Figure  10,  top,  we  show  a  different  initial  condition 
inside  O.  The  inpainted  image  (bottom)  is  practically  identical  to 
the  one  in  Figure  9,  thereby  showing  the  robustness  of  the  algo¬ 
rithm. 

Figure  1 1  shows  a  vandalized  image  and  its  restoration,  followed 
by  an  example  where  overimposed  text  is  removed  from  the  image. 
These  are  typical  examples  where  texture  synthesis  algorithms  as 
those  described  in  the  introduction  can  not  be  used,  since  the  num¬ 
ber  of  different  regions  to  be  filled  is  very  large. 


Figure  7:  Synthetic  example:  Q  is  shown  in  white.  Topology  is 
not  an  issue,  and  the  recovered  contours  smoothly  continue  the 
isophotes. 


Figure  12  shows  the  progressive  nature  of  the  algorithm:  several 
intermediate  steps  of  the  inpainting  procedure  are  shown  for  a  detail 
of  figure  11. 

Finally,  Figure  13  shows  an  entertainment  application  The 
bungee  cord  and  the  knot  tying  the  man’s  legs  have  been  removed. 
Given  the  size  of  9  a  2-level  multiresolution  scheme  was  used. 
Here  it  becomes  apparent  that  it  is  the  user  who  has  to  supply  the 
algorithm  with  the  masking  image,  since  the  choice  of  the  region  to 
inpaint  is  completely  subjective. 

5  Conclusions  and  future  work 

In  this  paper  we  have  introduced  a  novel  algorithm  for  image  in¬ 
painting  that  attempts  to  replicate  the  basic  techniques  used  by  pro¬ 
fessional  restorators.  The  basic  idea  is  to  smoothly  propagate  infor¬ 
mation  from  the  surrounding  areas  in  the  isophotes  direction.  The 
user  needs  only  to  provide  the  region  to  be  inpainted,  the  rest  is 
automatically  performed  by  the  algorithm  in  a  few  minutes.  The 
inpainted  images  are  sharp  and  without  color  artifacts.  The  exam¬ 
ples  shown  suggest  a  wide  range  of  applications  like  restoration 
of  old  photographs  and  damaged  film,  removal  of  superimposed 
text,  and  removal  of  objects.  The  results  can  either  be  adopted  as 
a  final  restoration  or  be  used  to  provide  an  initial  point  for  manual 
restoration,  thereby  reducing  the  total  restoration  time  by  orders  of 
magnitude. 

One  of  the  main  problems  with  our  technique  is  the  reproduction 
of  large  textured  regions,  as  can  be  seen  in  Figure  14.  The  algorithm 
here  proposed  is  currently  being  tested  in  conjunction  with  texture 
synthesis  ideas  to  address  this  issue. 

The  inpainting  algorithm  here  presented  has  been  clearly  mo¬ 
tivated  by  and  has  borrowed  from  the  intensive  work  on  the  use 
of  Partial  Differential  Equations  (PDE’s)  in  image  processing  and 
computer  vision.  When  “blindly”  letting  the  grid  go  to  zero,  the  in¬ 
painting  technique  in  equations  (5)-(ll),  naively  resembles  a  third 
order  equation  (not  resulting  from  a  variational  gradient  descent 
flow,  see  below).  Although  adding  regularization  terms  might  make 
this  high  order  equation  stable  and  allow  a  formal  analysis,5  to  the 
best  of  our  knowledge  the  complete  understanding  of  such  type  of 
equations  is  not  just  beyond  the  scope  of  this  paper  but  beyond  the 
current  state  of  mathematical  knowledge  (although  results  for  other 
high  order  equations,  which  might  be  relevant  for  image  processing 
as  well,  are  available,  e.g.,  [22]).  In  addition,  the  discrete  algorithm 
here  introduced  uses  as  “boundary  conditions”  both  gray  values  and 
isophote  directions  present  at  the  narrow  band  surrounding  the  re¬ 
gion  to  be  inpainted.  This  is  done,  as  explained  above,  performing 
the  computations  in  the  dilated  region  .  All  this  important  infor¬ 
mation  needs  to  be  incorporated  to  such  an  equation  as  “extra”  (and 

5  Note  that  numerical  implementations  of  PDE’s  intrinsically  regularize 
the  equations. 


fundamental)  boundary  conditions.  Nevertheless,  this  suggests  the 
investigation  of  the  use  of  lower,  second  order,  PDE’s  to  address  the 
inpainting  problem.  Moreover,  having  these  equations  to  be  gradi¬ 
ent  descent  flows  of  variational  formulations  permits  the  inclusion 
of  the  narrow  band  information  simply  by  changing  the  limits  of 
integration.  One  possible  approach  to  follow  is  to  first  reconstruct 
the  isophote  directions  0  (a  planar  vector  field),  and  from  them  the 
corresponding  gray  values  I.  Once  the  isophote  directions  are  re¬ 
constructed,  the  gradient  direction  is  given,  and  we  can  reconstruct 
the  gray  levels  by  finding  an  image  consistent  with  these  gradi¬ 
ent  directions.  Although  this  can  be  done  using  a  one  dimensional 
transport  equation,  the  existent  theory  limits  the  possible  gradient 
fields,  thereby  disqualifying  edges.  We  can  take  then  a  variational 
approach  and  find  the  image  I  minimizing  J*  (|  VJ|  —  0  •  V/)c?0, 
where  0  is  the  reconstructed  gradient  direction  inside  Oe  (an  ad¬ 
ditional  term  that  penalizes  the  deviation  from  the  real  image  in 
the  narrow  band  Oe  —  O  can  be  added  as  well).  The  minimal  I 
is  obtained  via  gradient  descent  flows  (second  order  PDE).  There¬ 
fore,  what  is  left  is  to  reconstruct  the  vector  field  0.  This  can  again 
be  achieved  using  a  variational  formulation  on  0  with  the  proper 
boundary  conditions  to  guarantee  the  continuity  of  the  vector  field 
across  dfl  (see  also  [13,  21]).  Both  variational  formulations  can 
actually  be  coupled  to  simultaneously  recover  the  isophotes  (vector 
field)  and  gray  values  via  the  corresponding  gradient  descent  flow, 
obtaining  a  set  of  coupled  second  order  PDE’s.6  Preliminary  results 
in  this  direction  are  promising  and  will  be  further  investigated  and 
reported  elsewhere  [23].  One  of  the  advantages  of  this  approach  is 
that  formal  results  regarding  uniqueness  and  existence  of  the  solu¬ 
tions  can  be  shown,  thereby  reducing  the  necessity  to  rely  mainly 
on  experimental  validations  as  done  until  now  with  the  basic  image 
inpainting  algorithms  reported  in  the  literature. 
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Figure  9:  The  user  defines  the  region  to  inpaint  (here  shown  in  red). 


Figure  10:  From  a  different  initial  condition  inside  0  we  arrive  at  a  very  similar  result.. 


Figure  1 1 :  Restoration  of  a  color  image  and  removal  of  superimposed  text. 


Figure  12:  Progressive  nature  of  the  algorithm:  several  intermediate  steps  of  the  reconstruction  of  figure  1 1  (detail). 


Figure  13:  The  bungee  cord  and  the  knot  tying  the  man’s  feet  have  been  removed. 


Figure  14:  Limitations  of  the  algorithm:  texture  is  not  reproduced. 


